Monte Carlo evaluation of path integrals for the nuclear shell model 
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We present in detail a formulation of the shell model as a path integral and Monte Carlo tech- 
niques for its evaluation. The formulation, which linearizes the two-body interaction by an auxiliary 
field, is quite general, both in the form of the effective 'one-body' Hamiltonian and in the choice 
of ensemble. In particular, we derive formulas for the use of general (beyond monopole) pairing 
operators, as well as a novel extraction of the canonical (fixed-particle number) ensemble via an 
activity expansion. We discuss the advantages and disadvantages of the various formulations and 
ensembles and give several illustrative examples. We also discuss and illustrate calculation of the 
imaginary-time response function and the extraction, by maximum entropy methods, of the corre- 
sponding strength function. Finally, we discuss the "sign-problem" generic to fermion Monte Carlo 
calculations, and prove that a wide class of interactions are free of this limitation. 
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I. MOTIVATION AND ORGANIZATION 



Exact diagonalizations of shell model Hamiltonians in the Os-ld shell demonstrated [|| that the shell model can yield 
^ ' an accurate and consistent description for a wide range of nuclear properties in different nuclei, if the many-body basis 
is sufficiently large. However, the combinatorial scaling of the many-body space with the size of the single particle 
basis or the number of valence nucleons restricts such exact diagonalizations to small nuclei or nuclei with few valence 
particles M. 



In facing the general challenge of developing non-perturbative methods to describe strongly interacting many- 
body systems, various quantum Monte Carlo schemes have been proposed as an alternative to direct diagonalization 
-««^ Among them, the auxiliary-field Monte Carlo method is suitable for addressing interacting fermions ji). This 

method is a Monte Carlo evaluation of the path-integral obtained by the Hubbard-Stratonovich transformation ||] 
I , of the imaginary-time evolution operator. The many-body wavefunction is represented by a set of determinantal 
^ • wavefunctions evolving in fluctuating auxiliary fields. The method thus enforces the Pauli principle exactly, and 
^ , the storage and computation time scale gently with the single-particle basis or the number of particles. Auxiliary 
^ ■ field methods have been applied to condensed matter systems such as the Hubbard Model yielding important 

' information about electron correlations and magnetic properties. 

• i-H 

k> ^ In this paper, we discuss the application of auxiliary-field Monte Carlo techniques to the nuclear shell model. This 
; I involves a two-body interaction more general than the simple on-site repulsion of a Hubbard model. Our goal is 
to develop new methods for extending the applicability of the shell model, as well as to investigate the powers and 
limitations of auxiliary-field Monte Carlo methods for general fermion systems. 

We have previously published a letter with selected results for static observables in sd- and /p-shell nuclei. This 
paper serves to give the details of the implementation, introduce a method for calculating dynamical correlations 
and strength functions, demonstrate the method with simple (sd-shell) nuclei, and discuss several important issues 
that arise in the implementation. We also explore the limitations imposed by the negative contributions in the 
path-integral, referred to as the "sign problem" in the literature. 

Our presentation is organized as follows. We begin in Section II by using the Hubbard-Stratonovich (HS) trans- 
formation to write the imaginary time evolution operator exp(— as a path integral. This requires that the 
Hamiltonian H be cast as a quadratic form using an appropriate set of operators. We discuss in Section HI two ways 
in which this can be accomplished, using either particle density or pairing operators. The imaginary time evolution 
operator can be used to extract information about the system at finite temperature or in its ground state (/3 ^ oo). 
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The formulae for obtaining static obscrvables are presented in Section IV where methods for handling the canonical 
ensemble are also introduced. We then discuss in Section V the extraction of the strength functions for operators 
from the imaginary time response function. In Section VI, we briefly describe the computational algorithms for im- 
plementing our methods and present selected results of calculations for sd-shell nuclei. In the final section, we address 
the sign problem and also discuss a class of non-trivial interactions that give rise to a positive definite path-integral 
for some nuclei. 



II. THE IMAGINARY TIME EVOLUTION OPERATOR 

Given some many-body Hamiltonian H, we seek a tractable expression for the imaginary time evolution operator: 

U = exp{-pH) . (2.1) 



Here, /3 has units of inverse energy and can be interpreted as an imaginary time. (Here and throughout, we take 
h = 1 and measure all energies in units of MeV.) It is also clear that U can be interpreted as the partition operator 
for temperature /3~^. We will refer to U as the evolution operator hereafter. The operator H is usually a generalized 
Hamiltonian and might contain terms beyond the true Hamiltonian, such as — /xTV in the grand-canonical ensemble 
or —LoJz if we are 'cranking' the system. 

There are two formalisms for extracting information from the evolution operator: the "thermal" formalism (on 
which we will concentrate) and the "zero-temperature" formalism (to which the thermal formalism reduces in the 
limit /3 oo). In the thermal formalism, we begin with the partition function 



Z = Trexp(-/3iI) , 



(2.2) 



and then construct the thermal observable of an operator O: 



0) = -Tr 



(2.3) 



Here, the trace Tr is over many-body states of fixed (canonical) or all (grand-canonical) particle number. In the 
zero-temperature formalism we begin with a trial wavefunction ipo and use the evolution operator to project out the 
ground state, assuming that tpo is not orthogonal to the ground state. The expectation value of O is then given by 



O 



lim 





cxp(- 


|i?)Oexp( 


-IH) 




{^0 


eM-PH) 





(2.4) 



In this section, we describe how to write [/ in a form that allows Eqs. (2.3) or (2.4) to be evaluated. 



A. Path integral formulation of the evolution operator 



We restrict ourselves to generalized Hamiltonians that contain at most two-body terms. The Hamiltonian H can 
then be written as a quadratic form in some set of 'convenient' operators Oa- 

H = Y,^»da + lY.^adl, (2.5) 

a a 

where we've assumed that the quadratic term is diagonal in the Oa- The meaning of 'convenient' will become clear 
shortly, but typically it refers to one-'body' operators, either one-particle ('density') or one-quasiparticle ('pairing'). 
The strength of the two-body interaction is characterized by the real numbers V^. 
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For H in the quadratic form (2.5), one can write the evolution operator [/ as a path integral. The exponential is 
first split into Nt 'time' slices, (3 = NtAf3, so that 



U= cxp(-A/3i/) 



Nt 



(2.6) 



Then we perform the Hubbard-Stratonovich (HS) transformation on the two-body term for the n'th time slice to give 



exp 



(-A/3h)~ y°° ]Jda„„ (^M^^ ^ exp|-A/3|^^i|y„|aL+eaOa + SaK.a„„a„^| 



(2.7) 



where the phase factor Sq is ±1 if < and is ±i if Va > 0. Each real variable (Tan is the auxiliary field associated 
with Oa at time slice n. 

The approximation (2.7) is valid through order A/3, since the corrections are commutator terms of order (A/5)^. 
The evolution operator is then 



U 



exp 



{-A(3h) 



P^* [a]G{a) exp {-A(3h, (r^J) ... exp (-A(3h, (n)) 



(2., 



where the integration measure is 



Nt 



n— 1 ot 



the Gaussian factor is 



and the one-body hamiltonian is 



G(a)=exp(-^i|yjaL) 

\ an / 



(2.9a) 



(2.9b) 



(2.10) 



It is sometimes convenient to employ a continuum notation. 



U 



J V[a]expl-^jyTY,\Va\<yliT 



Texp 



dr hair) 



where T denotes time-ordering and 



V [a] = lim X>^' [a\ 

Nt^oo 



Texp 



/ dTh„{T)\ ^ lim TT exp (-A/3^^(r„ 
7o / ^ 



(2.11) 



(2.12) 



(2.13) 



In the limit of an infinite number of time slices Eq. (2.8) is exact. In practice one has a finite number of time 
slices and the approximation is valid only to order A/3. The case of only one time slice is known as the Static Path 
Approximation (SPA); previous work on the SPA and its extensions can be found in Refs. § and §. 
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Rewriting the evolution operator as a path integral can make the model space tractable. Consider the case where 
the Oa are density operators. Then Eq. (2.1) is an exponential of two-body operators; it acts on a Slater-determinant 
to produce a sum of many Slater-determinants. In contrast, the path-integral formulation (2.8) contains only ex- 
ponentials of one-body operators which, by Thouless' theorem [|lO[, takes a Slater-determinant to another single 
Slater-determinant. Therefore, instead of having to keep track of a very large number of determinants (often many 
thousands for modern matrix-diagonalization shell model codes such as OXBASH |l^), we need deal only with one 
Slater-determinant at a time. Of course, the price to be paid is the evaluation of a high-dimensional integral. However, 
the number of auxiliary fields grows only quadratically with the size of the single particle basis while the corresponding 
number of Slater-determinants grows exponentially. Furthermore, the integral can be evaluated stochastically, making 
the problem ideal for parallel computation. 



B. Monte Carlo evaluation of the path integral 



Formulating the evolution operator as a path integral over auxiliary fields reduces the problem to quadrature. For 
a limited number of auxiliary fields, such as in the SPA with only a quadrupole-quadrupole interaction, the integral 
can be evaluated by direct numerical quadrature. However, for more general cases (typically hundreds of fields), the 
integral must be evaluated stochastically using Monte Carlo techniques. 



Using the one-body evolution operator defined by 
we can write Eq. (2.3) or (2.4) as 

'o) = 



JV[a] G(a)(0(a))c(a) 
JV[a] G(a)C(a) 



For the zero-temperature formalism 



C(a) = (Vo Ua{P,0) 



and 



while for the thermal formalism (canonical and grand-canonical), 

C(a)EEtr[C/„(/3,0)] , 

and 



Tr 



O) = 



(2.14) 



(2.15) 



(2.16) 



(2.17) 



TrC/,(/3,0) 



(2.18) 



(2.19) 



To evaluate the path integral via Monte Carlo techniques, we must choose a normalizable positive-definite weight 
function W^, and generate an ensemble of statistically independent fields {ffi} such that the probability density to 
find a field with values ai is Wo-. . Defining the 'action' by 
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(2.20) 
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the required observable is then simply 



/DM (6) e-^. *E,(6).i>. 



where is the number of samples, 



$i = e-^VVFi (2.22) 



and Si = S„. , etc. Ideally W should approximate exp(— 5) closely. However, exp(— 5) is generally not positive and 
can even be complex. In some cases, $i may oscillate violently, giving rise to a very small denominator in Eq. (2.21) 
to be cancelled by a very small numerator. While this cancellation is exact analytically, it is only approximate in the 
Monte Carlo evaluation so that this 'sign problem' leads to large variances in the evaluation of the observable. 

There are several possible schemes for both the choice of W and the sampling of the fields. We typically choose 
W = |exp(— <S)| and generate the samples via random walk (Metropolis) methods. 



III. DECOMPOSITIONS OF THE HAMILTONIAN 



^ To realize the HS transformation, the two-body parts of H must be cast as a quadratic form in one-body operators 
Oa- As these latter can be either density operators or pair creation and annihilation operators (or both), there is 
considerable freedom in doing so. In the simplest example, let us consider an individual interaction term, 

H = a\a2a4as , (3.1) 

where aj, at are anti-commuting fermion creation and annihilation operators. In the pairing decomposition, we write 
(using the upper and lower bracket to indicate the grouping) 

H = ajaj 0403 (3.2a) 
= ^(aittj +0304)^ - ^(a|a| - 0304)^ + ^[ajaj, 0304] . (3.2b) 

The commutator is a one-body operator that can be put directly in the one-body Hamiltonian ha- The remaining two 
quadratic forms in pair-creation and -annihilation operators can be coupled to auxiliary fields in the HS transformation. 

In the density decomposition, there are two ways to proceed: we can group (1,3) and (2,4) to get 

H = 0^03 0204 —aj 04(523 (3.3a) 

= -a|a4i523 + ^[a|a3,a2a4] -|- ^(a|o3 -|- ala^f - ^{a\a3 - 0304)^ , (3.3b) 

or group (1,4) and (2,3) to get 

H = — ala4 a^as +a\a3624 ■ (3.4a) 

= ala3^24 - ^[a{a4,a2«3] - ^{a\a4 + alaa)'^ + ^{a\a4 - 0^03)^ • (3.4b) 

Again the commutator terms are one-body operators, but now the quadratic forms arc squares of density operators 
that conserve particle number. We refer to Eq. (3.3) as the 'direct' decomposition and Eq. (3.4) as the 'exchange' 
decomposition. 
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For any general two-body Hamiltonian, we can choose the pairing or density decompositions for different parts of 
the two-body interaction. Moreover, even within a pure density break-up decomposition, there is still freedom to 
choose between the direct and exchange formulations. Although the exact path integral result is independent of the 
scheme used, different schemes will lead to different results under certain approximations (e.g., mean field or SPA). 
The choice of decomposition will also affect the rate of convergence of our numerical result as iVt ^ cxd, as well as the 
statistical precision of the Monte Carlo evaluation. Most significantly, it affects the fluctuation of $ in Eq. (2.21) and 
thus determines the stability of the Monte Carlo calculation. (See Section VI below.) 

In the application of these methods to the nuclear shell model, it is particularly convenient to use quadratic forms 
of operators that respect rotational invariance, isospin symmetry, and the shell structure of the system. We introduce 
these in the following subsections for both the density and pairing decompositions. 



A. Density decomposition 

We begin by ignoring explicit isospin labels and by writing the rotationally invariant two-body Hamiltonian as 

^2 = ^ E E ^^(«^' C'^) E AlMicib)AjMicd) 
abed J M 

= i E E + + ^-'^f' K7^(°^' E AM{'^b)A.jM{cd) (3.5) 

abed J M 

where the sum is taken over all proton and neutron single-particle orbits (denoted by a, 5, c, d) and the pair creation 
and annihilation operators are given by 

4M(a&) = E i3amajbm,\JM)al^^a]^^^^ = -[a]^ x a]f^' (3.6a) 
AjM{ab) = ^ {jamajbmb\JM)aj^,n^aj^„i^ = [aj^ x a^J"'*^ . (3.6b) 

rriaTiib 

The Vj{ab,cd) are the angular-momentum coupled two-body matrix elements of a scalar potential V{fi,f2) defined 
as 

Vjiab,cd) = ([^,Jfi) X ^jAr2)f^'\Virur2m,Ari) x ^,.(^2)]^*^) , (3.7) 
(independent of M) while the anti-symmetrized two-body matrix elements Vf{ab, cd) are given by 

Vfiab, cd) = [(1 + Sab){l + S,d)r'/^ [Vj{ab, cd) - {-l)^^+^'-'Vj{ab, dc)] . (3.8) 

Before continuing discussion of the density decomposition, we note that the two-body Hamiltonian for fermion 
systems is completely specified by the set of anti-symmetrized two-body matrix elements Vf{ab,cd) that are the 
input to many standard shell model codes such as OXBASH jll]. Indeed, we can add to the V_f{ab,cd) any set of 
(unphysical) symmetric two-body matrix elements Vf {ab,cd) satisfying 

vfiab, cd) = (-l)^'=+^''-^F/(a6, dc) , (3.9) 

without altering the action of H2 on any many- fermion wave function. However, note that although the Vf {ab,cd) 
do not alter the eigenstates and eigenvalues of the full Hamiltonian, they can (and do) affect the character of the 
decomposition of H2 into density operators, as is shown below. In what follows, we define the set of two-body matrix 
elements Vj^ (ab, cd) that may possess no definite symmetries as 

Vf{ab, cd) = Vf{ab, cd) + Vf{ab, cd) , (3.10) 
allowing us to write the two-body Hamiltonian as 
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^2 = ^ E E + + Scd)f' Vfiab, cd) E A\^{ab)AjM{cd) . (3.11) 



4 

abed J M 



To obtain the density decomposition of H2, we perform a Pandya transformation to recouple (a, c) and (&, d) into 
density operators with definite multipolarity 

pKM{ab) = E {jamajbmb\KM)al^^^aj,mt , (3-12) 

where Sj^ma = (~l)''"^'""'*joTOa- Then H2 can be rewritten as 

H2=H'2 + H[\ (3.13) 

^2 = ^ E E^^^"^'^*^) PK-M{ac)pKM{hd) , (3.14) 

a6cc/ if M 

where the particle-hole matrix elements of the interaction are 



£;K(ac,6d) = (-l)^''+^=E(-l)'(2-^+l){j-^ t i 



and -ff( is a one-body operator given by 



x^Vj'{ab,cd)^{l + Sab){l + Scd) , (3.15) 



-^1 = E ^ad/5oo(a, d) , (3.16a) 

ad 



with 



'ai = -i E E(-l)''^'"^'' (2-^ + l) ^2j! + l ^-^''^"^' ''^)V(l + M(l + M • (3.16b) 

Note that adding symmetric matrix elements is equivalent to using the exchange decomposition for some parts of the 

interaction. The freedom in choosing the combinations of direct and exchange decomposition is then embodied in the 
arbitrary symmetric part of the matrix elements Vj^ . 

Introducing the shorthand notation i = {ac),j = (bd), we can write Eq. (3.14) as 

^2 = iEE^^(^'J')(-l)'''5xM(*)/5K-M(j) . (3.17) 

ij K 



H2 as 



Upon diagonalizing the matrix EK{i,j) to obtain eigenvalues Xxa and associated eigenvectors VKa, we can represent 

^2 = J E >^K{a){-l)''pKM{a)pK-M{a) , (3.18) 



2 Ka 



where 

PKm{oi) = ^PKM{i)VKa{i) ■ (3.19) 
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Finally, if we define 



QKMia) = (pKMia) + i-l)^ PK-M{a)) , (3.20a) 

1/2(1 +dMo) 



PKM{a) = (pKMia) - (-l)^pK-M(a)) , (3.20b) 

V2(l + Omo) 



then becomes 



Ka M>0 

This completes the representation of the two-body interaction as a diagonal quadratic form in density operators. We 
then couple auxiliary fields (Tkm{o:) to Qkm and TKM{ct) to Pkm in the HS transformation. (The latter are not to 

be confused with the "imaginary time" r.) 

In the treatment thus far, protons and neutrons were not distinguished from each other. Although the orig- 
inal Hamiltonian H2 conserves proton and neutron numbers, we ultimately might deal with one-body operators 
PKMicLpT^n) and PKAiiinjbp) {n,p subscripts denoting neutron and proton) that individually do not do so. The 
one-body hamiltonian hn appearing in the HS transformation then mixes neutrons and protons. The single-particle 
wavefunctions in a Slater determinant then contain both neutron and proton components and neutron and proton 
numbers are not conserved separately in each Monte Carlo sample; rather the conservation is enforced only statistically. 

It is, of course, possible to recouple so that only density operators separately conserving neutron and proton 
numbers {pKM{ap,bp) and pKM{an,bn)) are present. To do so, we write the two-body Hamiltonian in a manifestly 
isospin-invariant form, 

^2 = ^ E E [(1 + ^ab){l + Scd)f^ VfT{ab,cd)Al^.^^^{ab)AjT;MT.icd) , (3.22) 

abed J 

where, similar to the previous definition (3.6), the pair operator is 

4T;MTMb)= E Oama,i6m6|JM)(it„,iib|TT,)at^„^,^at^„^,^ . (3.23) 

Here {^,ta), etc. are the isospin indices with ta = — ^ for proton states and ta = ^ for neutron states, and (TTz) 
are the coupled isospin quantum numbers. The two-body Hamiltonian can now be written solely in terms of density 
operators that conserve the proton and neutron numbers. Namely, 



where 



with 



and 



H2=H[+H'2, (3.24) 



ad t=p,n 



= -7 E E(-1)^^'"^'^(2J + 1) \ yj^^^i(ab, bdy{l + 5a,){l + 5,a) , (3.26) 
4V J V2ja + 1 



^2 = JE E EKT{ac,bd)[pK,T{i) X pk,tU)V^° ■ (3-27) 



2 

abed K,T=0,1 
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Here, we define pkm,t as 

PKM,T = PK M,p + {—^)^pKM,n > (3.28) 

and the Ek,t are given by 

£;K,T=o(ac,6rf) = (-iy^+>^(-l)-^(2J+l)|j ^|x/(TT«TTM 



1 

X - 
2 



yfr=i(a6, cd) + ^(yj^T=o(«&, erf) - yfr=iK erf)) 



(3.29) 



i?i^,T=i(ac,6rf) = -(-l)^^+^=X^(-l)'^(2J + l)|j ^1 V(l + M(l + M 

X J (^Olr=o(a^ erf) - y^T=i(«?', erf)) . (3.30) 

In tliis isospin formalism, since AjT-MTz{(ib) = {—iy''~^^''~'^~^^AjT;MTz{ba), the definitions of the symmetric and 
antisymmetric parts of Vj^{ab,cd), Vjrp{ab,cd), and V_^{ab,cd) become 

Vj^'^iab, cd) = i [Vf^iab, cd) ± (-l) '+^»+^'^+^-iyj^(&a, erf)] . (3.31) 

Note that these expressions allow less freedom in manipulating the decomposition since we have to couple proton 
with proton and neutron with neutron in forming the density operators. Also note that EK,T=o{ac, bd) — EK,T=i {ac, bd) 
is an invariant related only to the physical part of the interactions, {Vfrp^^ + V^/V=o)- We can choose all Ek,t=i 
to be zero in the above (by setting Vfrp^^ = V^j^^q) leaving Ek,t=o completely determined by the physical matrix 
elements. In that case, we can halve the number of fields to be integrated. However, while introducing the isovector 
densities requires more fields, it also gives more freedom in choosing the unphysical matrix elements to optimize the 
calculation. 

If we now diagonalize the EKrihj) as before and form the operators 

QkmM = -^=1=={pkM,t{(^) + {-l)"^ PK-M,T{a)) , (3.32) 
-^2(1 + dM,0) 

PKM,T{a) = ' A PKM,T{a) - {-l)'^pK-M,T{a)) , (3.33) 

^2(1 + dM,o) 

the two-body part of the Hamiltonian can finally be written as 

^2 = ^EE^^^('^) E {Qkm,t{^)+Pkm,t{^)) ■ (3.34) 

KT a M>0 

In this decomposition, the one-body hamiltonian ha- of the HS transformation does not mix protons and neutrons. 

We can then represent the proton and neutron wavefunctions by separate determinants, and the number of neutrons 
and protons will be conserved rigorously during each Monte Carlo sample. For general interactions, even if we 
choose nonzero Ek,t=i matrix elements, the number of fields involved is half that for the neutron-proton mixing 
decomposition, and the matrix dimension is also halved. These two factors combine to speed up the computation 
significantly. In this sense, an isospin formalism is more favorable, although at the cost of limiting the degrees of 
freedom embodied in the symmetric matrix elements Vf . 
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B. Pairing decomposition 



In nuclei where the pairing interaction is important, it is natural to cast at least part of the two-body interaction as 
a quadratic form in pair creation and annihilation operators. We demonstrate this for the case where the Hamiltonian 
is written in the isospin formalism. Upon diagonalizing V^{ab,cd) in Eq. (3.22), we can write 



^2 = V 



J2 AjT(a) ^It-mt. {a)AjT-MT, (a) , (3.35) 

JTa MT^ 

where 

^It-.mt^o:) = X!^''^«*^*)^^'T,MT,(«) • (3-36) 

i 

Separating A^A into commutator and anticommutator terms, we have 

H2 = H^ + H[; (3.37) 

^'^ = IY1 ^■^^(") E [^JT;MT, (a), ^JT;MT, (a)] ; (3.38) 



JTa MT^ 



H2 = IY1 ^•^^(") 12 {AjT;MTAc^)\^JT;MTA(^)} ■ (3.39) 



JTa MT. 



Clearly, H[ is a one-body operator that can be combined with Hi. The remaining two-body term can be written as 
a sum of squares by defining 

QjT;MT, {a) = [^JT;MT, (") + ^JT;MT^ {oi)j (3.40a) 

PjT;MT, (a) = --^ {At;MT, (") " ^JT;MT, («)) , (3.40b) 

SO that 

^2 = ^ E E (Q't;MT. («) + P!t;MT. (")) • (3-41) 
JT MT^ 

As in the density decomposition, wo can then couple the a and r fields to Q and P, respectively. 

Note that in the pairing decomposition, the one-body Hamiltonian /i(t) used in the path integral is a generalized 
one-body operator that includes density, pair-creation and pair-annihilation operators. The wavefunction is then 
propagated as a Hartree-Fock-Bogoliubov state, rather than as a simple Slater determinant. 

In this decomposition, neutrons and protons are inevitably mixed together in the one-body Hamiltonian hg- (consider 
the Q,P terms for T = 0). In fact, h^r also does not conserve the total number of nucleons; rather the conservation is 
only statistical after a large number of Monte Carlo samples. 

For simplicity, we have described how to decompose the Hamiltonian solely in density operators or solely in pair 
operators. However, it is straightforward to mix the two decompositions with the choice depending on the type of 
interactions involved. Consider the 'Pairing plus Quadrupole' model, namely 

H2 = -gP^P -\xQ-Q, (3.42) 
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where , P are the monopole pair creation and annihilation operator and Q is the quadrupole- moment operator, 

^^-E«"^- (3.43) 

a 

0M-E("I^Ml/3)4a/?- (3.44) 

Naively, it would be most convenient to use a pairing decomposition for the pairing interaction and a density decompo- 
sition for the quadrupole interaction. This would require only 8 fields for each time slice in the HS transformation. The 
pure pairing or pure density decompositions would be much more complicated, as the former would require rewriting 
the quadrupole interactions in terms of pair operators and the latter would require rewriting the pairing interaction 
in terms of density operators. The numerical examples given in Section VI below will illustrate the behavior of such 
Hamiltonians under different decompositions. 



IV. CALCULATION OF STATIC OBSERVABLES 



In the previous sections, we expressed the evolution operator for a two-body Hamiltonian as a path integral over 
auxiliary fields in which the action involves only density and pair operators. In this and the following section, we 
show how to extract observables from this path integral. The merit of this method will be clear from the compact 
formulae involved, which require handling only relatively small matrices of dimension or 2Ns, depending on 
the type of decomposition used. We derive formulae for three different approaches that use the evolution operator 
to obtain information about the system: the zero-temperature formalism, the grand-canonical ensemble, and the 
canonical ensemble. The zero-temperature and grand-canonical ensemble methods have been applied, with the density 
decomposition, to other physical systems such as the Hubbard Model |^. However, we believe that the canonical 
ensemble treatment presented here is novel. 

The zero-temperature approach can be used only to extract ground-state information. On the other hand, the 
grand-canonical ensemble allows finite-temperature calculations, but fluctuations in particle number can be very 
significant in a system with a small number of nucleons. Thus, the canonical ensemble is particularly important in 
nuclear systems, where the particle number is small and shell structure is prominent. The grand-canonical ensemble 
yields information averaged over neighboring nuclei, which can have very different properties. 

The canonical ensemble is more difficult than the other two approaches in two respects. First, the canonical (fixed- 
number) trace of Ua- is more difficult to compute than the wavefunction overlap of the zero-temperature formalism or 
the grand-canonical trace. Second, observables are more difficult to extract since Wick's theorem does not apply. We 
suggest two different methods to handle the canonical ensemble: the activity expansion and the integration over real 
coherent states. The coherent-state integration method can be applied to calculate canonical trace of any particle 
number involving only a negligible increase in computation time relative to the zero-temperature and grand-canonical 
approaches. Unfortunately, its utility is limited by the sign problem. On the other hand, the activity expansion, 
which is well-suited (and numerically stable) for calculating nuclei with a small number of valence particles or holes, 
is less susceptible to the sign problem. (Yet a third method for extracting the canonical trace, based on Fourier 
analysis, shows good "sign" statistics and is numerically stable for mid-shell calculations; it is therefore suitable when 



the activity expansion fails. Details of this method will be published elsewhere |12|.) 



For each approach, we also derive the general formulae when pair operators (as well as density operators) are present 
in the single-particle Hamiltonian; i.e., when a pairing decomposition is used for some or all of the interaction. As far 
as we know, there is no known general formulae for the pairing decomposition formalism. Using the fermion coherent 
state formalism [ p^ , we derive in Appendix A a set of formulae for the calculation of observables that are similar to 
the well-known formulae for a pure density decomposition. Thus, our methods can be extended to calculations using 
general one-body operators in the HS transformation by only doubling the dimension of the matrices involved. 
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A. Zero-temperature formalism 



We begin with 'zero-temperature' observables. The trial wavefunction ■00 in Eq. (2.4) is, in principle, any state not 
orthogonal to the ground state. In practice, it is most conveniently a Slater-determinant. If we have Ng m-scheme 
orbitals available and iV„ indistinguishable particles, ipo is constructed from Ny single-particle wavefunctions, each of 
which is a vector with Ns components, and we write V'o as a Slater-determinant of a x Ng matrix, 'S'q. 

In the pure density decomposition, consider the one-body Hamiltonian h = Mija\aj (sums over the indices i, j from 
1 to Ng are implicit). The evolved wavefunction, 

|0(A/3)) = exp (-AM) IV'o) (4.1) 

is then a Thouless transformation of the original determinant, the new state being represented by the matrix 
exp(— A/3Af )^'o- We can therefore represent the product of evolution operators by Ns x Ng matrices 

exp (-A/3/i„) . . . exp (-A/3^i) ^ exp (-A/3Af„) . . . exp (-A/3Mi) . (4.2) 



Let us now consider the overlap function in Eq. (2.16). Let Ua{T2,Ti) be the matrix representing the evolution 
operator Ua{T2, ti). Choosing some value t of the imaginary time at which to insert the operator O, we introduce the 
right and left wavefunctions iPr.l{t) defined by 



IV'fl(r)) -t/.(r,0)|Vo) 



(4.3a) 
(4.3b) 



Then the required overlap in Eq. (2.16) is 



det *[*_R, , 



(4.4) 



where 



vI>«(t) = U,{t, 0)*o ; *l(t) = Ulifi, t)*c 



(4.5) 



are the matrices representing V'i?, and ip^. Note that if there are two distinguishable species of nuclcons — protons 
and neutrons — and we use a decomposition that conserves the numbers of neutrons and protons, then there is a 
separate determinant for each set of single-particle wavefunctions and the total overlap is the product of the proton 
and neutron overlaps. 

With the basic overlap in hand, we now turn to the expectation value of an operator O for a given field configuration 
(i.e., Eqs. (2.15) and (2.19)). By Wick's theorem, the expectation value of any A^-body operator can be expressed as 
the sum of products of expectation values of one-body operators. Hence, the basic quantity required is 



alaa ) = 



*1 



(4.6) 



(A straightforward derivation of this expansion can be found in [p^.) Again, when the decomposition separately 
conserves proton and neutron numbers, the expectation values for proton and neutron operators are given by separate 
matrices. 

These formulae can be extended to the case where pair operators are present in the one-body Hamiltonian; i.e., 
when h is of the form 



(4.7) 
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where Q, A, and A are general Ns x Ng complex matrices. Our first task is to find a simple expression for 
{''Po\Ua{P,0)\4'o) where Ua{P,0) = Ilnii 6xp(— A/J/i^^)- Using the fermion coherent state representation and Grass- 
man algebra, we derive the expressions in Appendix A. Here we simply state the results. 



If the trial wavefunction ipo is a quasi-particle vacuum, such that /3i\4>o) = where 0i = ^ ■ UijUj + VijO.^-, then 



(Vol n exp(-A/3/i„)|Vo) = det 



n=l 



{v* n*)t/<.(/3,0)( ;t 



(4.8) 



where 

U„{0, 0) = exp(-A/3M;vJ • • • exp(-A/3Mi) , (4.9) 
is the matrix representing the evolution operator C/o-(/3, 0), and M„ is the 2Ns x 2Ns matrix representing 

'^- = Ia*'aI "'-ef)^ ("0) 



Here, the evolution operator U„{t2,ti) is represented by a 2Ns x 2Ns matrix and the many-body wavefunction is 
represented by a 2Ns x matrix, independent of the number of particles present. In analogy to Eqs. (4.3-4.6), we 
can write 



f>a(/3,0) Vo) = (VL|V'fl) = det 



^-p(-T^ETr[e„] 



(4.11) 



where 



^ J7.(r, 0) ( J ) , *^ ^ [/t (/3, r) ( J ) 



(4.12) 



To calculate the expectation value of a generalized one-body operator, we proceed as in the pure density case. Let 

aj=Oj, i = l,...,Ns (4.13a) 



1, 



(4.13b) 



Then the one-quasiparticle density matrix is 

{i>L \Pab\ tpR.) 



- [*H(*i*ii)-'*i 



1 






.ba~ 2 







ba 



ab 



■ Sab 



(4.14) 



The final stop follows from the properties of U in Eq. (A20). Note that both the overlap and the Green's function 
are similar to those of the density decomposition. However, the enlarged dimension of the representation matrices 
causes the overlap to be the square-root of a determinant rather than just a determinant. We know of no simple 
way to determine the sign of the square root. Computationally, it can be traced by watching the c^'olution of 
{'fpo\Ua{T,0)\tjjo) as Ua is built up from to /3 (Appendix B), although at the expense of increased computation. 
Also, the linear dimension of the matrix used in the calculation is twice that for the pure density decomposition case 
so that the pairing decomposition is more computationally demanding. However, when the interaction has a strong 
pairing character, it has the potential for a more effective Monte Carlo sampling, and offers greater freedom in the 
decomposition, which might be used to mitigate the sign problem. 
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B. Grand-Canonical Ensemble 



For the grand-canonical ensemble, the trace in Eq. (2.18) is a sum over all possible many-body states with all possible 
nucleon numbers, and a chemical potential in H is required. For the pure density decomposition, the many-body 
trace is given by 



C(a) = Tr f/,(/3, 0) - det (1 + i7,(/3, 0)) 



(4.15) 



which can be proved by expanding the determinant [|l4| . 

For the expectation value of a one-body operator, one can recapitulate exactly the argument of the zero-temperature 
development and obtain 



ba 



(4.16) 



We have extended these formulae to decompositions that involve pairing operators (Appendix A). The results are 
given in terms of the 2Ns x 2Ns matrices Mn, U{/3,0) representing the Hamiltonians ft.„ and the evolution operator 
U{P,0), (Eqs. (4.9,10)) namely 



Tr 



C/(/3,0) =det[l-f t/(/3,0)]5exp 



A/3 



Tr[e„: 



(4.17) 



To motivate this formula, consider the simplest case of one time slice where U = exp(— /i) and h is hermitian and 
in the form (4.7). Then, h can be diagonalized by a HFB transformation 



(at a) 



e A- A^ 

A - -e^ 



(a at) + -Tr[e] 



(4.18a) 
(4.18b) 



where 



u V 
V u* 



(4.19) 



and 



e 








L-A^ 






) 


/ei 


... 





...\ 


£2 


... 










... -ei 










... 


-£2 










■J 



(4.20) 



In the diagonal form, Tr[exp(— /i)] can be identified easily as 



)ei^']e-*T^[®l = W v/(l + e"'0(l + e^Oe"* 



lTr[e] 



det[(l 



e-^)ile-iTr[e] 



(4.21) 



due to the invariance of the determinant with respect to similarity transformations. Here, as in the overlap formulae 
for zero-temperature approach, the grand-canonical trace is given as the square root of a determinant, so that the 
evolution of the sign is important (Appendix B). The formula for the observables can be calculated from 
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(ata,)^ =-[(1 + [/(/?, 0))-^C/(/3,0) 



1 

1 



{bboxl + U(/3,0)y^U(/3,0) 



1 

1 



ab 



(l + f/(/3,0))-^J7(/3,0) 



6a 



(4.22) 



C. Activity expansion for the canonical ensemble 



As mentioned in the beginning of this section, the grand-canonical ensemble may lead to large fluctuations in the 
particle number for systems with few particles, and so is particularly ill-suited for small nuclear systems, and although 

the particle number does not fluctuate in the zero-tcmpcraturc approach, that formalism can only give ground-state 
results. The canonical ensemble is therefore important for studying thermal behavior, as well as the ground state of 
large systems. 

In the canonical ensemble, we have to find the trace of Ua{(3,0) over all states with a fixed particle number 
(actually fixed proton and neutron numbers). We discuss two methods of achieving this: the activity expansion 
presented here and the integration over real coherent states presented in the following subsection. 

Consider first the case when only density operators are present in the one-body Hamiltonian h. Prom the grand- 
canonical ensemble formulae (4.15), one can sec that the canonical trace is just the s\im of all the Ny X Ny sub- 
determinants. More explicitly, we consider the activity expansion: for some parameter A, let 

Z{l3, A) = tr A*" e-^" = ^ A^- Zn, (/3) . (4.23) 

In our matrix representation, 

tr A*;7<,(/3,0) = det(l + XU) . (4.24) 

If Eq. (4.24) is expanded in powers of A, the coefficient of A^" is just the canonical trace of UaiP, 0) over Ny particles, 
so that det(l -|- XU) is the generating function for the canonical trace. Thus, 

ZnM = J P[a]G(a)Cjv„(a) , (4.25) 

where 

det (1 + Af/,(/3, 0)) = J2 >^''^CnM) ■ (4.26) 



The trick now is to find simpler expressions for C.^^ , instead of doing the explicit sum over the determinants. To 
do this, write 

det(l + XU) = exp Tr ln(l + XU) = exp ( ^ A" Tr [[/"] J . (4.27) 

This expansion converges to the generating function because Z{f3, A) is a polynomial in A of finite order (i.e., Ny can 
be at most Ng). The coefficient of A" in the exponential is readily found. For a given particle number Ny, we need 
only calculate matrix traces up to Tr[l7^''] and the coefficient of A^" can be extracted accordingly. 

For one-body expectation values, using again the grand-canonical trace as the generating function for calculating 
observables and collecting all terms with coefficient A^" , we arrive at 
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(ala,)^^^^ = J2i-ir-' {U-\a CN-nia)/CN{a) . (4.28) 



71=1 



The expectation value of two-body operators, {pabpcd{<y)) n^-, is nontrivial as Wick's theorem must be modified, but 
again one can simply collect the terms with coefiicient A^" and obtain 

{a\a,alaa)^^^^ = E E [(-1)"+" {ULUZ - U^JJZ) C^.-™-„(a)/C^, (a)] 

n=l \m=l 

+^6c(-l)'^''-'C/3aCiv„-„(<7)/Cjv„(a)| . (4.29) 

This form of the activity expansion works best for < Ns/2 (mostly empty model spaces). When Ny > Ns/2 
(mostly filled spaces), it is more efficient to expand in the activity of the hole states. In this case, we define 

Z^{p, A) = tr[A^=-*e-^*] , (4.30) 

and as in Eq. (4.24), 

Tr[A^=-*C/(/3, 0; cr)] = det[A + U] . (4.31) 
The coefficient of A^ is the partition function for TV holes, 

det[A + U]= det[U] exp(^ -(-l)"-U"Tr[l7-"]) = ^ X^Cn , (4-32) 



n 

n>l N 



the expectation value of a one-body operator is 



N 



{Pij)N holes = E(-l)"^7i"^^^f-n/^^^f ' ^^.33) 



n=0 



and the expectation of a two-body operator is 

N 

{PijPkl)N holes = E ^ifc(~^)"'^(i"CiV-n/Civ) 



n=0 

N 



-(^ (-i)'"+"(f/-,"f/rr - uT,"u-ncN-m-n/CN ■ (4.34) 



m=0 



When pairing operators are involved, we use Eq. (4.17) for the grand-canonical trace, which becomes the generating 
function for the corresponding canonical trace. As an illustration, we display the formula for the expansion in particle 
activity, 



Tr 



A^i7 



det 



Al 








u 



A 2 



det 



Al 











■0 




522 J 



A^ 



= det 



512- 
522 



det 



+ A 



512 
522 



5" 
521 



det 



512- 
522 



exp 



iTr 
2 



In 







1 / 1 \n 

= det [5^2] = exp - [F"] (-l)"-i — 



+ XY 



(4.35) 
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where the definitions of S^^, S^^, S^^, S"^^ are obvious, and 

The expectation values of one- and two-body operators can then be derived as in the pure density decomposition. 

D. Canonical ensemble via coherent states 

We make use of the operator 



ph 



which can be shown [|5| to be a resolution of unity in the Hilbert space of iV^, fermions occupying levels. Here, 
Xph are the real integration variables, \X) are the real coherent states 

|X) =exp(^Xp,,ata,,)|0), = 1, . . . , 7V„; p = A^, -f 1, . . . , iV, , (4.38) 

hp 

|0) is the A'^-fermion state with the levels 1, . . . , occupied and C is a normalization constant, 

Yl Cp, c, = .-^=/2r(| + i)/r(^ + i) . (4.39) 

p=N^ + l ^ ^ 

The canonical trace can then be cast in a form of expectation with integration over the variables Xpi^. For any 
operator O, 

tr[6] = C [ dX ^^'^'^^ , (4.40) 

J det{l+XTX)^+^ 

where dX = Y[ph ^-^ph ■ The thermal canonical expectation is then 



jV[a]G{a)dX{X\U,{p.,Q)d\X)/dei{l+X^Xr^l^+^ ^^^^^^ 



J V[a]G{a)dX{X\U„{/3, 0)|X)/ dct(l + XTX)N^/^+^ 

We can compute the overlap and the observable in (4.41) as in the zero-temperature formalism, but we now have 
to do the Monte Carlo integration over the fields Xph as well. However, the number of these fields scales only as A^^, 
which is about the same as the number of a fields in one time slice, so that introducing the coherent-state integration 
is not much more computationally expensive than the zero-temperature formalism. The advantage of coherent states 
is that we do not have to find Tr[L'"(/3, 0)"]. However, the integration of the coherent states may aggravate the sign 
problem, as will be discussed in Sections VI and VH. 



V. DYNAMICAL CORRELATIONS 

In the previous section, we discussed how to extract the expectation value of a one-body operator (O) and, by use of 
Wick's theorem and its extension to the canonical case, equal-time two-body operators {AB), etc. A great deal of the 
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interesting physics, however, is contained in the dynamical response function (O^(t)O(O)) where 0{t) = e ^HtQ^tHt^ 
In our calculations, we evaluate the imaginary-time response function (O^(ir)O(O)) and from it deduce the associated 
strength function, f^{E). 

In the zero temperature formalism, the strength function is 

f^iE) ^Y.^{E-Ef+ Emm? , (5.1) 

/ 

where i is the ground state, while in the canonical or grand-canonical ensemble, 

foiE) ^\Y.^{E-Ef+ E,)e-^^^ \ (/|0|z) p . (5.2) 

Thus, the imaginary-time response function is related to the strength function by a Laplace transform, 

1*00 

i?(T)^(Ot(,r)d(0))= / f^{E)e-^^dE. (5.3) 



—00 



Recovering the strength function from R by inversion of the Laplace transform is an ill-posed numerical problem. 
Different methods have been proposed to surmount this difficulty . We resort to making the best "guess" for the 

strength function via the Classic Maximum Entropy (MaxEnt) method ||l^,|l^, which was first introduced to recover 
strength functions in Monte Carlo calculations by Silver et al. |17]| and has since been widely used in similar condensed 
matter simulations. It is in essence a least-squares fit biased by a measure of the phase space of the strength function. 

In MaxEnt methods, the function to be minimized is ~ where is the usual square of the residuals, 
S is the entropy of the phase space (not to be confused with the action in the auxiliary field path integral), and 
a is a Lagrange multiplier. In Classic MaxEnt a is determined self-consistcntly. The method, described briefly in 
Appendix C, also can yield error estimates for the extracted strength function. 



To calculate the imaginary time response, consider i^A{iT)B{Q) ) and write the thermal expectation as a path 
integral: 

Trfe^^'^^^^^Ae^^-^Bl 
(A(zt)B(O)) = - ^ ^ 



Tr[e-/3^^] 



ltr[&„(/3,T)Ai/„(r,0)B] 



(5.4) 



/l?[a]G(a)Tr[C/,(/3,0)] 
Upon defining 

OAt) = Ua{T, oy^OU^iT, 0) (5.5) 
we have an expression suitable for Monte Carlo evaluation, 

{A(iT)Bm = „ . Tr[[/.(ff,0)] ,g g . 



/l?[a]e-^(-)(i.(r)g.(0)) 



(5.6b) 



We now proceed to find Oct(t). For the purpose of illustration, we show the formulae for the pure density decom- 
position (formulae for the general case can be derived similarly) . For the simplest case when O = a| or O = , it can 
be shown that 10 



18 



4^nm = ^[C/.(nA/3,0)-i]5-a] (5.7a) 

j 

aUnAp) - J2 UainAp, 0)^0^ . (5.7b) 

3 

In this way, the creation and annihilation operators can be 'propagated' back to r = 0. Any operator Oa-ir) can be 
first expressed in a„{T) and a^o-(''") and therefore can be propagated back and expressed in a and a^. For example, 
suppose O = Cija\aj is a one-body operator. Then 

da{T) - Ua{T, Q)-^dU„{T, 0) (5.8) 

- [(7,(r,0)-iC(7<,(T,0)]yaJa, . (5.9) 
Thus, the response function can be measured in the same way as the static observables. 



VI. COMPUTATIONAL DETAILS AND ILLUSTRATIONS 



A. Monte Carlo Methods 



Monte Carlo evaluation of the path integral requires a weight function. We have tried two different choices for the 
weight function, each with advantages and disadvantages. 

One choice for the weight function is a Gaussian, with the the static mean-field solution as the centroid, and the 
widths given by the RPA frequencies. Thus, much of the known physics is embodied in the weight and the Monte 
Carlo evaluates corrections to the mean-field -I- RPA approximation. Further, it is simple to efficiently generate 
random field configurations with a Gaussian distribution. 

Unfortunately, Gaussian sampling has several disadvantages. First, finding the RPA frequencies can be extremely 

time consuming since we have to calculate and diagonalize the curvature matrix ^ — at the mean- field 

solution. (Here, a, 7 run through the number of quadratic terms in the Hamiltonian and i,j run through the number 
of time slices.) For any general interaction, the curvature matrix has a large dimension, N'^Nf. Second, the Gaussian 
has to be modified when there is spontaneous symmetry breaking in the mean fields (such as quadrupole deformation) . 
Otherwise, the Goldstone modes in the the RPA frequencies (e.g., zero eigenvalues corresponding to shape rotations) 
will destroy the normalizability of the weight function. Finally, multiple mean-field solutions well separated from each 
other can also pose a problem, so that multiple Gaussians may be needed. 

A more satisfactory route is to choose | exp(— 5)| — G{a)\({Ua{P, 0))| as the weight function and to use a stochastic 
random walk such as the Metropolis algorithm to generate the fields. The expectation of an observable is then given 
by Eq. (2.21) where 

= cxpham(5,)] . (6.1) 

The Metropolis algorithm is free of the disadvantages for the Gaussian weight function, in that it will eventually 
sample the entire space where | exp(— 5)1 is significant. The main disadvantages of Metropolis are its inefficiency as 
currently implemented and the "sign problem." Let us define the denominator of Eq. (2.21) by ($): 

<'^)^^E^^pHM'50]- (6.2) 

i 

If ($) ^ 1, the large fluctuations defeat any attempt at a Monte Carlo evaluation. This 'sign problem' will be 
addressed in detail in the examples illustrated below. 

For the Metropolis algorithm, we take random steps in the fields time-slice by time-slice, following a sweeping 
procedure introduced by Sugiyama and Koonin For the Monte Carlo results to be valid, one requires that the 
points in the random walk be both distributed as the weight function and be statistically independent. The first 
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requirement translates into starting the fields in a region of statistically significant weight; if the initial configuration 
has a small weight, a number of initial "thermalization" sweeps are usually needed to relax the fields to this region. 
The second requirement means that the walker must sweep through the fields many times to de-correlate the samples. 

The number of thermalization sweeps and dccorrclation sweeps increases with system size, and the choice of random 
walk procedure greatly affects the sampling efficiency. In the early stage of investigation, we allowed all fields Ga 
at one time-slice to change with equal probability within a certain step size. The acceptance is then determined 
by the ratio of the weight |exp(— 5)1 of the old and new configurations, and we adjusted the step size to give an 
average acceptance probability of 0.5. The calculations of srf-shell nuclei described below then needed approximately 
2000 thermalization sweeps and up to 200 decorrelation sweeps between samples. In our later work, we used another 
random walk algorithm, in which a fixed number of fields at one time slice arc randomly chosen for updating; those 
chosen are generated from the Gaussian distribution in Eq. (2.9b) while the others are kept at their previous values, 
and the acceptance is determined by the ratio of Q in the new and old configurations. This alternative algorithm 
needs only some 200 thermalization sweeps and 10 decorrelation sweeps; i.e., it is 10 times more efficient than the 
previous method. In addition, boundaries where | exp(— iS)| = can confine the first random walk algorithm, while 
they do not affect the second. 

The random walk can thermalize faster if it starts from a configuration where the weight W{a) is large. Usually 
we start from the static mean fields, given by 

croc,n = CTa, n= 1, Nt ; (6.3a) 

CTa = -SoiSlgn(ya){da)a ■ (6.3b) 



One can show that for canonical and grand-canonical calculations, the static mean field Oa is a saddle point of the 
weight function, 

(6.4) 

For the zero-temperature approach, we also choose the self-consistent field solution ct^, although Eq. (6.4) is not 
rigorously satisfied. This is preferable to starting the fields at zero, which may be far from configurations of significant 
weight. 

The mean field solution (6.3) can be found iteratively 

= -SaS\gn{Va){6a)ai ■ (6.5) 

For the systems we have tested, Eq. (6.5) usually converges within 100 interactions starting from CTo- = 0. For larger 
systems and at lower temperature, convergence is slow and sometimes unstable or barely stable. In that case taking 

CTa.i = -SaSign(Va)(6„)^=0 (6.6) 

for a starting configuration also leads to a shorter thermalization time than the choice of a = 0. 



B. Examples 



We now describe several examples to demonstrate our methods for calculating nuclear properties. Two major 
considerations arise in the implementation: the choice of decomposition scheme and the choice of ensemble. Different 
decomposition schemes involve different dimensions of matrices and numbers of fields, which control the computational 
speed. Also, the rates of convergence as A/3 are different and determine the number of time slices to be used. 
Finally, and perhaps most importantly, this choice also affects the "sign problem" associated with the statistical 
stability of the calculation. Different decomposition schemes will be compared in Examples 1 and 2 below. 

The choice of ensemble among zero-tempcrature, canonical, and grand canonical ensembles usually do not affect the 
issues noted above. Instead, it depends on the kind of properties to be calculated. The zero-temperature formalism 
with a good trial wavefunction is effective in projecting out the ground state and is suitable for calculating ground-state 
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static observables. For calculating finite temperature properties, the canonical ensemble is physically most relevant but 
also is most difficult to implement. Examples 1, 2, and 3 below demonstrate the grand canonical, zero-temperature, 
and canonical ensembles, respectively. 

Finally, we choose a particular nucleus, ^"Ne, to demonstrate the calculation of the dynamical responses for different 
operators and recover the strength function by the MaxEnt Method. The examples shown below were done with 3000 
to 6000 samples generated on the passively parallel Touchstone Gamma and Delta computers at Caltech. 



Example 1: Monopole pairing interaction in the sd- shell 



We have described the considerable flexibility in writing the two-body interaction in quadratic form; e.g., density 
versus pairing decomposition and direct versus exchange decomposition. To illustrate this flexibility, we consider 
protons only in the sd-shell [Ng = 12), and keep only the J = terms in Eq. (3.5); the values of Vj^q and single 
particle energies are taken from the Wildenthal interaction 0. We first recouple the Hamiltonian into a quadratic 
form in the density operators in Eq. (3.21); because all possible density operators are required, there are 144 fields for 
each time slice. We then use the pairing decomposition in Eq. (3.41); after diagonalization only 6 fields are required. 
Finally, we write the two-body interaction as H2 — ^-^2 + and decompose the first half using densities and the 
second half using pairs; the total number of fields in this case is 150. It turns out that this system is 99% free of the 
sign problem (i.e., Real(<i>i) > 99% of the time), independent of the decomposition. 

All three calculations were performed in the grand-canonical ensemble using a Gaussian weight function around 
the static mean field, at an inverse temperature of /? = 1 (here and henceforth, we measure all physical energies in 
MeV) and fixed chemical potential. The expectation value of the proton number, energy, and are given in Fig. 1 
as a function of A/3. As the number of time slices increases (i.e.. A/? 0), all three decompositions converge to the 
exact answer, demonstrating their mutual equivalence in the continuum limit. Note, however, the different rates of 
convergence. 



Example 2: Mg with schematic forces 



Next we consider sd-shell nuclei with a schematic Pairing -t- Multipole density interaction, where the multipole force 
is separable; it is the same interaction used in H: 



Tj = -l,0,l 

-^X2^P2,MP2 -m(-I)^'"^ - ^X4^P4.MP4,-Af(-1)^^ (6.7) 



M M 



where 



Pt, = ^ {^ti'^t2\'^T^)ajmtiajmt2 (6.8) 

jniti t2 



and 



PK,M = UK{jlj2)pK,M,T=o{jl,j2 



(6.9) 



This Hamiltonian was also decomposed in several ways. We first decomposed the pairing interaction as pair operators 
and the multipole force as density operators. In this way, the number of fields is kept to a minimum (only 21 per time 
slice). The pair-operator Ptj=o mixes protons and neutrons and therefore one matrix representing the mixed neutron 
and proton wavefunction is needed. In addition, the pairing decomposition requires matrices whose dimension is 2Ns, 
so that the matrices involved in this calculation are 48 x 48. 

We calculated ^"^Mg (4 valence protons and 4 valence neutrons) in the zero temperature formalism; i.e., using the 
evolution operator at large /3 to project out the ground state from a trial state ipQ. Since h is hermitian (here p 
has the property p]^ ^ ~ Pk,~m{^^)^^), in the static path C = {ipo\ exp{— f3h)\ipo) is always positive definite and 
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$i = 1. However, for calculations with more than one time slice ($) becomes very small, so that we can only obtain 
sensible results in the SPA. These results turned out to be extremely good, relaxing to the right energy and angular 
momentum (Fig. 2). However, the success of the static path is case-specific and not well understood. For example, 
we have also tried the case of just multipole interactions among 4 protons in the sd-shell (^°Mg), and find that the 
SPA relaxes to an energy 2 MeV higher than the ground state. 

In a second scheme, we transformed the pairing part of the interaction (6.7) into a quadratic form in the density 
operators. In this transformation, we used only density operators that conserve proton and neutron numbers (3.24- 
3.34), and chose all Ek.t=i elements in Eq. (3.27) to be zero so that only isoscalar density operators were present in 
the quadratic form. 

In this case, the interaction is in a much more complicated form due to the Pandya transformation of the pairing 
interaction, and one needs 144 fields per time slice (as compared to 21 fields needed in the first decomposition). 
However, an advantage lies in the separation of the Slater-determinants for protons and neutrons since only density 
operators that separately conserve neutron and proton numbers are present. In the sd-shell, the dimension of matrix 
involved is only 12 x 12. For this particular interaction, an even more desirable property of the second scheme is that 
the eigenvalues A^.a in Eq. (3.34) found by diagonalizing Ek in Eq. (3.27) satisfy sign(Aif ) = (—1)^+^. We prove in 
the next section that this property guarantees Ci'^) to be positive definite for even-even nuclei if a suitable trial state 
is chosen. This allows calculations with any number of time slices that are free from any sign problem. 

We performed the calculation using the zero-temperature formalism at different /3 and A/3 values, choosing first 
the Hartree Slater-determinant as the trial wavefunction. The SPA calculation and that with A/3 = 0.125 are shown 
in Fig. 2. We have also performed calculations at A/3 = 0.5 and A/3 = 0.25. These results are not shown but lie 
between the SPA and the A/3 = 0.125 results, and show a convergence to the result at A/3 = 0.125. At A/3 = 0.125, 
the energy relaxes to the right energy, whereas the SPA also relaxes, but to a much higher energy. 

We repeated these last calculations with a different trial wavefunction: a Slater determinant of the orbitals (j, m) = 
(|,±i), (|,±|). Although a different relaxation curve is traced out by the results at A/3 — 0.125, it also converges 
to the same exact result (Fig. 2). The choice of the trial state is therefore important for determining the rate of 
relaxation of the zero-temperature approach, although the final result is independent of the trial state as expected. 
In this case, although the Hartree state is lower in energy than the maximal prolate state (compare (H) at /3 = 0), 
it contains some high angular momentum components (compare (J^) at /? = 0), so that it relaxes more slowly and 
reaches the ground state at a larger value of (3. 

Example 3: Canonical calculations of'^^Ne and '^'^Mg 

Next, we demonstrate the canonical ensemble for the same interaction (6.7) using the pure density decomposition 
as described in Example 2. We calculate ^°Ne because it is small enough to allow for an exact diagonalization to 
give all the states of iJ, since we are concerned with both ground-state and thermal properties. The canonical trace 
for this path- integral in also positive definite (see Section VII) , allowing the calculations to be done using many time 
slices. 

The results for calculations with A/3 = 0.25, 0.125, and 0.0625 are shown in Fig. 3. The convergence as a function 
of A/3 is also apparent. Note that particular attention should be given to the extrapolation at high temperature. 
However, it is not hard to increase the number of time slices to decrease the error at high temperature. For example, 
for /3 = 0.5, A/3 = 0.0625 amounts to only 8 time slices. Similar results on ^'^Mg in the canonical ensemble are shown 
in Fig. 4. The relaxation to the ground state can be compared with the zero-temperature result in Fig. 2; however, 
now the results at small values of (3 are also physically significant. 

In the calculation of ^°Ne, the activity expansion method is numerically stable. However, instabilities appear for 
s(i-nuclei when the number of proton or neutron valence particles (or holes) is greater than 4. The instability is 
signalled by the deviation of (Np) and (Nn) from integers. In those cases, the expansion in Eq. (4.27) or (4.35) gives 
the canonical trace as the small difference between very large numbers, so that mid-shell nuclei cannot be calculated 
directly by those equations. We have developed an alternative Fourier expansion technique to extract the canonical 
trace that gives satisfactory results . 

The real coherent-state method for the canonical trace for ^''Mg gives the same results as the activity expansion. 
However, is not always unity due to the need to integrate over the coherent states in Eq. (4.40) (as will be 
explained in the next section). It changes from 0.70 to 0.23 for /3 changing from 0.5 to 1.0 at A/3 = 0.25. 

We have also studied rotating nuclei in the canonical ensemble by adding a Lagrange multiplier to the Hamiltonian, 
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H' = H — LoJz, where Jz is the z-component of the angular momentum and lo is the cranking frequency. The 
calculations at /3 = 1 for ^°Ne with A/3 = 0.125, 0.0625, and 0.03125 are shown in Fig. 5, where the convergence to the 
exact results can be seen. The moments of inertia fitted from the three sets of data are 5.10, 5.30, and 4.95, compared 
to 4.74, the value from the exact curve. By adding the term linear in j^, we break the time reversal symmetry of h{a), 
which is related to the sign properties of the weight function C(c)- {^) decreases from 1 to 0.55 when to is increased 
from to 1.5 at A/3 = 0.125 while it decreases from 1 to 0.52 at A/3 = 0.03125. 

Example 4: Response and strength functions for ^°/Ve 

Finally, we demonstrate calculation of the imaginary-time response functions and the reconstruction of the strength 
functions. The calculation for ^°Ne is done by the activity expansion method in a pure density decomposition; the 
Hamiltonian is that of Eq. (6.7). The canonical ensemble is more suitable than the zero temperature formalism for 
measuring the dynamical response, because in the latter case many "boundary" time slices are needed to project out 
the ground state on both the left and the right, and extra time slices would have to be introduced in the middle to 
measure the response. In contrast, the cyclic property of the trace allows full use to be made of every time slice in the 
canonical ensemble. We choose /3 = 2.5 (the system has already reached the ground state at this low temperature, as 
can be seen from Fig. 3) and calculate the response functions at A/3 = 0.125 and A/3 — 0.0625 for several one-body 
operators: the isovector- and the isoscalar-quadrupolc operators Qv = Qp — Qn, Qs = Qp + Qn, and the isoscalar and 
isovector angular momentum operators Js = Jp + Jn, Jv = Jp — Jn- We choose this latter 1+ operator purely out of 
convenience; Js is just the total angular momentum, which we verified to produce a constant response equal to ( J^), 
which follows from the rotational invariance of H. 

The response functions are shown in a semi-log plot in Fig. 6(a,c,e). For these Hermitian operators, R{t) is 
symmetric about r = /3/2, so it is shown only for r < (3/2. The slope of the plot is approximately the energy of the 
dominant strength. The Jy and responses relax more rapidly than does the Qs response, indicating that the two 
isovector operators couple to states with higher excitation energy than does the Qs operator. Since ^°Ne has a J = 0, 
T = ground state, the states excited by an operator will carry the J and T quantum numbers of the operator. The 
plots are therefore consistent with the existence of a low-lying 2+ state. 

The MaxEnt reconstructions of the most probable strength function for the different one-body operators are shown 
in Fig. 6(b,d,c). The reconstruction is performed for each response hmction measured at each A/3 value. The figures 
show the convergence in A/3 of the resulting strength functions to the exact strength function. Note the movement 
of the peaks towards the exact position and also the decrease in the widths as A/3 decreases. 

The average n'th moments of the strength function, defined by 

Mr, = Y^oj^f(uJi); (6.10) 

i 

can be found in the Monte Carlo integration over all the possible distributions of the fi. Their uncertainties can be 
determined similarly. The 1®* moment (Mi) and the 2"'* moment (M2) of the strength functions are listed in Table 1 
for different operators and different A/3. The extrapolated moments (A/3 0) and the exact results for the ground 
state transitions are also shown. 

The single-particle pick-up and stripping response functions for different orbitals arc given in a semi-log plot in 
Fig. 7(a,c,e) and Fig. 8. The strength functions extracted from these responses are related to the excitation spectrum 
in the neighboring nucleus with one additional particle or hole. The stripping and pick-up responses are the same 
for protons and neutrons as the ground state of ^"Ne is isoscalar, and the final states have the angular momentum of 
the single nucleon that is added or removed. We see from Fig. 7 that both the pick-up and stripping responses for 
the j = I orbital and the pick-up response for the ) = ^ orbital converge to the exact results as A/3 becomes small; 
the MaxEnt reconstruction of the corresponding strength functions in Fig. 7(b,d,e) also show a convergence to the 
exact results. The extracted moments are listed in Table 1. However, the responses for the j — | orbital show an 
anomalous behavior: they arc close to the exact result at r = 0, and then, with a sudden change in slope, follow the 
responses of the j = | orbital. A possible reason is that angular momentum is not conserved sample-by-sample in the 
calculation, but rather only statistically. The -/ = | states for ^^Ne and ^"Ne nuclei are much lower in energy than 
the J = \ states (because the i = f orbital is lower than the J = | orbital by 5.6 MeV), so that if a small admixture 
of the J = I states "leaks" into the intermediate states for the ) = | response, it will dominate with an exponentially 
growing correlation function. (The J = ^ orbital is much closer to the J = f orbital in comparison, only 0.8 MeV 
higher.) 
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VII. THE SIGN PROBLEM AND FUTURE DEVELOPMENTS 



We have seen in the examples above that the sign properties for the function G{a)C{cr) are crucial to the stability 
of the calculation, and may frustrate attempts to apply the Monte Carlo path integral to any general two-body 
interaction. In general, if we choose an arbitrary two-body interaction and arbitrarily decompose it into pair or 
density operators, (<I>) vanishes rapidly as /? increases or as the number of time slices Nt increases at a fixed p. For 
example, with the Wildenthal interaction and for any calculation for (3 > 1 with more than three time slices, the 
noise due to the sign completely swamps the signal. This "sign-problem" is well-known in auxiliary-field path-integral 
calculations |20| and fermion quantum Monte Carlo more generally. If there is no explicit symmetry to enforce the 
positivity of Co-, ($) decays exponentially as a function of /3, and the problem is more severe for smaller A/3, so that 
it is very difficult to calculate low-temperature properties. 

Only a handful of interacting fermion systems are known to give rise to positive definite path-integral: the one- 
dimensional Hubbard Model, the half-filled Hubbard Model, and the attractive Hubbard Model at any dimension and 
filling 1^. We will show in the following that an important class of interactions for the nuclear shell model has a 
positive-definite path-integral representation for even-even nuclei. It arises from the symmetry between time-reversed 
single-particle orbitals and may serve as a starting point in understanding and controlling the 'sign problem'. 

We first define the 'time-reversed' partner of each creation and annihilation operator to be 

O'j^m — CLj.fh — ( 1)^ (^j,~m (7.1a) 

Note that 

due to the spin-half statistics. 

The class of Hamiltonians that give rise to a positive-definite path integral are of the form 

H = -]^'^XaPaPa+ eaPa+ e*^Pa (7.3) 
a 

where Xa > 0, Cq, can be generally complex, and pa is a general density operator of the form 

Note the requirement of a negative coupling constant for the density operator with its time-reversed partner and the 
time-reversal invariant form of the remaining one-body part. 

For application to the shell model, we refer to Eq. (3.17), so that in the density decomposition, 

PKAiia) - PK-M{a) . (7.5) 

The requirement of a negative coupling constant for the density operators leads to a sign rule for the Xkoj namely 

Sign(AKa) = (-1)^+' . (7.6) 



As we will show below, the Hamiltonian (7.3) gives rise to a one-body hamiltonian in the path integral, h„, that is 
symmetric in time-reversed orbitals. Time-reversed pairs of single-particle orbitals thus acquire complex-conjugate 
phases in the propagation, guaranteeing a positive definite overlap function C. 

After the HS transformation on Eq. (7.3), the linearized Hamiltonian is 

h = - ^Xa iio-a + iTa)Pa + {cTa " iTa)Pa) + BaPa + e*^Pa , (7.7) 
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so that p and p couple to complex-conjugate fields. (If some part of the interaction does not satisfy (7.6), then there 
will be terms in h that are of the form i{(7a + iTa)pa + *(o'a ~ iTa)Pa, and the above statement is not true.) 

We represent the single-particle wavefunction by a vector of the form 



jm 
jrri 



, m > , (7.8) 



with states with m > in the first half of the vector and their time reversed orbitals in the second half. Due to Eq. (7.3) 
and the fact that time reversed operators are coupled to complex conjugate fields, the matrix Mi representing hi has 
the structure 

^i={-B* A*) ' (^-^^ 
and one can easily verify that the total evolution matrix, 

U = l[ exp(-MiA/?) = ( _Q* p. ) , (7.10) 

is of the same form. Here, A, B, P, Q are matrices of dimension Ns/2. One can show that this matrix has pairs of 
c„.„p,» .,e„»,„s .... »,„«vee,e™c.„. («) .„d ( :X l - ^ w..e„ . ,e.l, i. . 



V I \ u 



doubly degenerate since the two eigenvectors are distinct. 

For the grand-canonical ensemble, the overlap function is given by 

r / M ^^^^ 

C = det l+f = n(l + e,)(l + e*)>0. (7.11) 

L \ / J j 

If only particle-type (neutron-proton) conserving density operators are present in Eq. (7.3), each type of nucleon is 
represented by a separate Slater-determinant having the structure (7.10), and therefore C = Cp ^ Cn > 0, since (p > 
and Cn > 0. 

In the zero-temperature formalism, if the trial wavefunction for an even number of particles is chosen to consist of 
time-reversed pairs of single-particle states, 

where a, b are matrices with dimension (^ x ^), then ^qU^q is a, Ny x matrix with the structure (7.10), and 
the overlap function again satisfies 

C = det[*Sf/*o] > . (7.13) 

If only particle-type conserving operators are present, then time-reversed pairs of trial wavefunctions can be chosen 
for both protons and neutrons in an even-even nucleus, giving rise to ^ = x Cn > 0. Note that while the overlap 
function is positive definite in the grand-canonical ensemble for any chemical potential (and thus any average number 
of particles), it is true only for an even (or even-even if involving a type-conserving decomposition) system in the 
zero-temperature formalism with a suitable trial wavefunction. 

In the canonical ensemble for N particles, a fixed-number trace is involved and therefore 

(7.14) 

\ i / il^i2^...iK 

We have found empirically that ^ is positive definite for even-even systems, although we lack a rigorous proof. A 
special case of the Hamiltonian (7.2) exists in which the overlap function is positive definite also for odd-odd N = Z 
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nuclei. The required condition is that only isoscalar-dcnsity operators arc present in Eq. (7.2). This leads to a further 
symmetry that protons and neutrons couple to the same field in Eq. (7.7) and therefore the evolution matrices satisfy 
Up = Un = U. In the zero-tcmpcratmo formalism, if we choose the trial wavefunction for proton and neutrons to be 
time reversed partners of each other, so that 




(7.15) 



then 

Cp = det[*S_pC7*o,p] 

= det[a^P a + a^Q b - b^Q*a + b^P*b] = Q , (7-16) 

and ( = Cp Cn > 0. On the other hand, for the canonical ensemble, one can prove that Eq. (7.13) is real. Given 
that Cn = Cp) C is a square and is therefore positive. 

For a system with general form (7.3), if we perform the canonical trace by integration over real coherent states, 
then there is an extra operator, exp{Xphaj,ah), multiplying the evolution operator. Time-reversed states couple to 
different real fields in the extra operator, and a sign problem arises, as seen in Example 3 above. 

Cranking also causes the sign to depart from unity. In cranking, the Lagrange multiplier term wJ^ is added to h in 
Eq. (7.7), destroying the property that time-reversed operators are coupled to complex-conjugate numbers (because 
Jz = —Jz)- Notice, however, that cranking with an imaginary Lagrange multiplier iujJz will preserve the time-reversal 

symmetry and will give rise to positive path-integral. 

In summary, for a Hamiltonian of the form (7.3), the above proof guarantees the overlap function to be positive 
for any nucleus in the grand canonical ensemble, and for even-even nuclei in either the canonical ensemble or zero- 
temperature formalism (with suitable trial wavefunction). It also guarantees positivity for odd-odd N = Z systems 
when only isoscalar density operators arc involved. 

The Hamiltonian (6.7) satisfies (7.6) upon decomposition of the paring interaction using density operators and 
involves only isoscalar operators, so that C, is positive for even-even and N ^ Z nuclei. For other systems, ($) 
decreases as a function of /?. At A/3 = 0.0625, ($) = 0.4 at /3 = 1.5 for ^^Na and ($) = 0.2 at /3 = 2.0 for ^SNa. Thus, 
even for odd-A nuclei, the sign properties of (6.7) are still much better than that of a general interaction violating the 
criteria (7.6). For example, using the Wildenthal interaction, ($) drops to several percent at /? = 1 for any nucleus. 

Arbitrary two-body interactions do not satisfy the sign rule (7.6), and when the rule is violated, ($) rapidly decreases 
as a function of /3. Note that monopole pairing plays an important role here. The pairing interaction can be written 
as 

ij 

It produces a constant shift in every Xa in Eq. (7.3), as may be seen also from the multipole decomposition of the 
pairing force: 

-gY,{-^fpKM{a)pK-Mia){-l)^ . (7.18) 

a 

Therefore, if the pairing interaction is strong enough compared to the remaining interactions, the sign rule can be 
satisfied. 



VIII. CONCLUSION 

We have developed a general framework for carrying out auxiliary-field Monte Carlo calculations of the nuclear shell 
model. In this framework we evaluate ground state or thermal observables, using pairing or density fields or both. 
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Although the use of pairing fields naturally embodies important aspects of the residual interaction, these calculations 
are more difficult due to the larger matrix dimension needed and the extra effort to keep track of the sign of the overlap 
function as the wavefunction is propagated. Furthermore, for calculations with multiple time slices, the Monte Carlo 
method with a pairing decomposition suffers from severe sign problems. However, pairing fields are suitable for carrying 
out static path or two-time-slice calculations where the linearized Hamiltonian is hermitian, thereby enforcing the 
positivity of the overlap function. This can be easily verified by observing that for hermitian h, ha, and hb, with real 
eigenvalues Ei,Ei^ and Ei^, 

tr[e-'5''] = ^6-''^^ >0, (8.1) 

i 

and 

trfe-^'^oe-^'^'] =^e-^^-(ia|e-^'"'|ia) >0. (8.2) 



In these cases, there is no sign problem and also there is no need to keep track of the evolution of the sign. 

For the density decomposition, we have found a class of interactions which give rise to a positive definite integrand 
upon the HS transformation. For these interactions, stable calculations can be carried out for many time slices to 
extrapolate to the exact results (A/3 — + 0). This class of interactions includes the phenomenological pairing-plus- 
multipole interaction used widely. We have carried out calculations with such interactions in the sc?-shell, demon- 
strating the power of the method in calculating both static and dynamical properties in the ground state and at finite 
temperature; high spin nuclei were also studied by cranking. The calculations converge to the exact results (as found 
by direct diagonalization) with increasing number of time slices. Although the nuclear wavefunction is not found 
explicitly in these calculations, many nuclear properties can be obtained. 

For general shell model interactions, it appears that the sign or phase property of the integrand is the major 
factor determining successful application of the Monte Carlo sampling. Successful calculations are usually confined to 
high temperature studies. We have demonstrated the freedom one has in the decomposition scheme of the two-body 
interaction and found that the behavior of the sign can be different in the various schemes. The next crucial step is 
to explore whether we can manipulate these degrees of freedom to enable stable calculations of nuclei using general 
forces. 
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APPENDIX A: DERIVATION OF THE OVERLAP FORMULAE FOR PAIRING FIELDS 



We consider operators of the form 

C/ = exp(-A/3/i7vJexp(-A/?/iArj_i) . . .exp(-A/3ft,i) , (Al) 

where each ht is a quadratic operator 

Ns 

ht= e(t)y a|aj + A{t),jalal + A{t),jaiaj . (A2) 

Without loss of generahty, we choose A, A to be anti-symmetric (A — — A-^, etc.). We follow the development of 
Berezin [Q, who considered the special case U — exp{—iht) with h Hermitian; we take the general case. 

We calculate the grand-canonical trace 

Ti-U = (A3) 

i 

^h e sum is over all states of all particle number) by using the fermion coherent state (FCS) representation of unity 

1 = / n^^"<^^p(-E^"^") 10(^1- (A4) 

Here are Grassman variables and the |^) are fermion coherent states. Then 

TtU = Jl[d(o.dCeM~J2i*c.iamm- (A5) 



We need the FCS representation of U. In what immediately follows we show that if U is the matrix representation 
of U, that is 

U = {jj2i JJ22 j (A6a) 
= exp(-M(A^f)A/3) exp(-M(iVt - l)A/3) . . . exp(-M(l)A/3), (A6b) 



where 



then 



_ f e{t) 2A{t) \ 



(^|[/|0=Cexp(l(C r)(^2! f2l)Q.]], (A8) 



with 



^11 ^ U22-Ijj21 ^ ^12 ^ _t722-l ^ (^g^) 
^21 ^ ([J22-1)T ^ ^22 ^ [7l2jj22-l ^ (^g^) 



and 
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Afi ^' 

C7 = det([/='=')iexp( — A^Tre(n)). (AlO) 

n=l 

In this case the trace becomes a Gaussian integral over Grassman variables; the result is given beginning with equation 
(A29) below. However, to come to that point we must derive equations (A8-A10). 

To this end, we employ the standard rules for operating on |^), (^| with ,a: 

m^uo = (Mia) 

maiio = (Mib) 

{^\aim = cm\o (auc) 

i^Wam = -^mo ■ (Aiid) 

Next, we derive expressions for UaU, a^jj, as linear combinations of t/aa, J/aJj. Then, with the ansatz (A8) for 

as a Gaussian in the Grassman variables ^, we use (All) to derive the elements B of the Gaussian given in (A9). 

To do this, we introduce the operators h, b (which are not necessarily Hermitian conjugates) 

6„ = tJ-^a^U, 6„ = U-^fltf/ . (A12a) 

Then 

ajj = Uba, alU = UK , (A12b) 

and we seek b, b as linear combinations of o^, a. 
Define 

a„(r) = e^^a„e-'*^ . (A13) 

Then 

^a„(r) = [/i,a„(r)] , (A14) 

and similarly for al^{T). Putting all the aa(r), al^{T) into a single vector, and using the representation (A2) for h, one 
finds 



(A15) 



with M given by (A7). 

Solving the differential equation (A15), 



and so in general 



«(;))=exp(-Mr)(« ) , (A16) 



= exp(AMi) . . . exp(A/3/ijvJ ( ) exp(-A/3ftjvJ • • • exp(-A/3/ii) 



= exp(-A/3Mi) . . .exp(-A/3MjvJ 

- I ^21 JJ22 I I 



(A17) 
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Then 

be = U^^a^a^ + U^^a-ya\ , (Al8a) 
be = U^^e^a^ + U^^e-ya\ , (A18b) 

where the summation on 7 is implicit. Upon inserting (AI8) in (Al2b), and using the ansatz (AS), one can straight- 
forwardly derive the B's in terms of the f/'s as given in (A9). 

Although we do not show it in detail, we note that B is antisymmetric {B^^^ — B^^, B^^^ = B^"^, etc.), which 
can be proved using 

lA ,^fo i\ 



and 



1 ; " 1 I = 



Now we must show the normalization C is of the form (AlO). To do so, we find a differential equation for C. Letting 

Un+i = eM-^l3hn+i)Un , (A21) 

we define 

Unit) = exp{-thn+i)Un (A22) 

so Un+1 = t^n(A/3). Taking the expectation value of (A22) between FCS's, invoking (All) and equating the parts 
independent of £,,£,*, one obtains 

^lnC„(t) = -Tr(A(n+l)B2\(i)) . (A23) 

Upon differentiating (AlO), one obtains 

I In Cn{t) = ^Tr^ In U^^,{t) - \ti-Q„+, , (A24) 

Using Un{t) — exp{—Mn+it)Un, one derives 

^U^\ ^ -2An+iU'\ + Ol+iU^^ (A25) 
at 

and (A24) becomes (A23). Thus Eq. (AlO) satisfies the differential equation for C, and it just remains to establish 
the overall normalization. This is found by choosing M = 0, so that [/ = 1, C = 1, and 

(^|J7|0 =exp(-^e:ea) (A26) 

a 

which is (Cl^). Thus we have established the form (A8) for (^|?7|^). 

The integral (A5) is straightforward (see Berezin ||l^; the magnitude is 

Tr?7 = C7det^fi^^ 522 ■ (A27) 
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The phase of TvU, though critical, is more difficult to obtain (see Appendix B for details). 

One can rewrite (A27) in a more compact form. The constant C contains the factor det(S^^)5, which can be 
written using 

det J7^2 = det ) = det ^''^^ ^ ) (-1)^^' . (A28) 

Upon introducing (A28) into (A27), performing some algebra, and using relationships from (A20), one arrives at 

AB ^' 

frU = det(l + U)^ exp( — ^Tr V @{n)) , (A29) 

where, again, U is the matrix in (A6) representing the evolution operator U. 

As for the density case, one can introduce an activity expansion to project out an exact particle number 

Tr(A^[/)=det(l+(^^ ? ^ t/)^ exp ^^Tre(n)j A^^/^ . (A30) 

Because the matrices 1 and U are of dimension 2Ns x 2Ns, one can write A^'/^ = det ^ ^ and (A30) becomes 

tr(A%)=det((j + xexp(^-^f:Tre(n)j. (A31) 

This can be expanded into a polynomial in A, which then gives the canonical ensemble. 

Finally, we give the expectation value of (V't|C/|V't). First we note that the vacuum expectation value (OlC/IO) is the 
term in (A31) independent of A, 

(0|[/|0)=det(j JJ,y^'exp||-^^Tre(n)j (A32a) 
= det((0 l)t/(5))%xp('-^^Tre(n)') . (A32b) 



Any quasi-particle excitations can be represented as Hartree-Fock-Bogoliubov vacua for properly dc^fincxl new quasi- 
particle operators, which corresponds to doing a similarity transformation on the matrices. If \ipt) is the vacuum to 
the quasi-particle annihilation operator i.e., 

AlV't) =0, A = UijUj + Vija], (A33) 

then 

(^W.)=det((0 J) (;)) xexp^-f^Trew) 

= detf(t,* u*)u(''l\]\^^(-^Y.^rQ{n)\. (A34) 
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APPENDIX B: SIGN OF THE OVERLAP FOR PAIRING FIELDS 



The formulae derived in Appendix A for calculating the overlap Ca{P) = {'4>l\U{P, 0; (T)\tpii) in the zero-tcmpcrature 
formalism and Tr[C/o-(/3, 0)] in the thermal formalisms all involve the square root of a determinant, leaving the phase 
of Cct(/3) undetermined by a factor of ±1. This ambiguity is irrelevant for the Monte Carlo random walk as we 
typically take \(^\ as the weight function, but the phase must be known unambiguously for calculation of observables, 
as important cancellations may result. 

We determine the phase by following the evolution of Ca{T) and its derivatives with respect to r as r goes from 
to (3. For example, if Cct(t) were purely real, zero-crossing (with C'(''') 7^ 0) would indicate a change in the phase by 
-1. The initial phase at r = is real and positive. Following the evolution is computationally expensive, but as most 
of the time is spent on the random walk, where the phase is irrelevant, the overall computational time is negligible. 

In what follows we give the formulae for up to fourth derivatives for each of the different formalisms. 

Greind Canonical Ensemble 

Define 

fit) = Ca(fcA/3 + t) = fT[e-f"'+''U4kAf3, 0)], (Bl) 

then 

k 

m = det[l + e-^-'+^^t/]^ exp(-^ ^Tr[e,] - ^lV[e,+i]) (B2) 

ln(/) = lTV[ln(l + e-^''+^*C/)] - ^fl ^^[0^1 - ^TV[e,+i] . (B3) 

Using the abbreviation M = Mk+i,& = ©fe+i, let 

G = (1 + e-^'U)-^e-^'U =!-(! + e-^*C/)-\ (B4) 

^ = _[1 + e-^*C/]-iMe-^*[/[l + e-^*C/]-i 

= -{l-G)MG. (B5) 



The derivatives of ln(/) can be expressed in terms of the matrices G and M, 

^ ^ -^Tr[MG] - iTr[e] (B6) 

92 ^ = ^Tr[M(l - G)MG] (B7) 

53 = = --Tr[MGM(l - G)M{1 - 2G)] (B8) 

54 = = -^TV[M(1 - G)MGM{1 - G)M{1 - 2G)] 

-iTr[MGM(l - G)MGM(1 - 2G)] - Tr[GM(l - G)M(1 - G)MGM] . (B9) 



Then Ca{kAp + t) is given by. 



Ca(fcA/? + t) = fit) = /(O) exp ( + + + + . . . ) . (BIO) 
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Zero temperature formalism 

In this case 



fit) = UkAp + t) = {i;t\e-'^''+HU{kAl3, 0; = (VL|e-^'=+^*|^ij) , (Bll) 

= det[*ie-^^+^**H]ie-^Et,^[0^1-|^[®^+il) , (B12) 

ln(/) = iTr[ln(*ie-^'=+^**R)] - ^fl ^^[9.] " In^k+i] ■ (B13) 



Let 



2 .=1 



G = e-^**R[*ie-^**H]-i*r , (B14) 
= GMG-MG = -(1-G)MG. (B15) 



Then 



= = -\Tr[MG] - IttIQ] , (B16) 

and so on and aU formulae are the same as in the grand canonical ensemble (B6-9), except that the matrix G is now 
different. In fact, in both cases G can be shown to be the matrix for the Green's function; i.e., 

Gij = (a]a,) (B17) 

where 

ai = ai,i = l,...,Ns (B18a) 

ai+N, =4,i = l,...,Ns . (B18b) 
Canonical Ensemble 

In Eq. (4.35), the undetermined sign involves only the vacuum expectation (0|C/(/?, 0)|0). Once it is determined, 
the sign for Tr[C/(/3, 0)] is known. We can use the equations for the zero temperature formalism to obtain the sign of 
(0|C/(/3,0)|0). 
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APPENDIX C: MAXIMUM ENTROPY EXTRACTION OF THE STRENGTH FUNCTION 



We use the Maximum entropy (MaxEnt) method to reconstruct the strength function from the response function. 
Here we give a brief description of the Classic MaxEnt Method; details can be found in the paper by GuU 

The MaxEnt method is a Bayesian approach for reconstruction of positive additive images / from noisy data. In 
our case the image is the strength function f{uj). The noisy data, Dj — dj +1]^, are the measurements of the response 
function R{t) at the discrete imaginary times, where dj = R{jA(3) and rjj is the noise in the data. In the absence of 
any data, the most probable image is chosen to be a default model to. Skilling |^8| proved that, in that case, the only 
consistent choice for the probability of an image / is determined up to a parameter a: 

pr{f) = exp(a5(/, TO))/Z(a, m) , (CI) 

where S is the entropy of image / relative to the default model. If the image is discretized to fi {i — 1, . . . ,r), then 

S{f, m) = Y.^f, - m, - /, ln(/,/TO,)) , (C2) 

3 

and 

/•oo /"OO 

Zsia,m)^ d''fl[f-iexip{aS{f))= d^u exp{aS{u^)) , (C3) 

Jo J ~oo 

where the last step follows from a change of variable Ui = ^/Ji- [Zs has no relation to the nuclear partition function.) 

In the presence of data, we gain some knowledge about the image. Assuming a Gaussian distribution of errors in 
the data Di, the probability for / is 

prif) = exp(a5 - ix'(/))^^ , (C4) 



where 



Gij = {r]ir]j) is the correlation matrix of the errors in the data, and 

Zl= Jd^'Dexpi-^x^f))- (C6) 

For a given choice of a, the most probable image can be found by maximizing aS — ^X^i giving rise to the term 
"Maximum Entropy" method. However, a is not predetermined. Rather, it is varied until the maximum of 

aS — is approximately equal to the total number of data Di. But this assignment of a is ad hoc and, according to 
Gull p9[ , usually leads to an underfitting of the data. In the classic MaxEnt, a is fixed by maximizing the probability 
of a given the data set Di and the default model to: 

pria\D,m)^ZQZs'Zz\ (C7) 

where 

Zq = d'f n exp(a^ - -x') = J d^u exp(a5 - -x") ■ (C8) 

After a is fixed, we can find the most probable / by maximizing aS — \x^ , or we can find the average / by Monte 
Carlo sampling using the integrand of (C4) as weight function. Information about the uncertainty in / can also be 
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obtained from Monte Carlo sampling. In the particular application at hand, we would like to know, for example, the 
uncertainty in the location of the peaks in /, or in the moments M„ ~ J f{Lu)uj^^cLj/ J f{uj)duj. 

Let us now return to the problem of extracting the strength function. Since the imaginary time response function 
R{t) is given at discrete times, we can only allow a limited amount of parameters in the strength function /, which 
we do by discretizing f{u>) to fi at oji. To allow for some smoothness, we choose a Gaussian function centered around 
each oJi, rather than a delta function. Given the imaginary time response function, we bound the range of lu over 
which the / is significant by tOmin and tOmax, and choose the Ui to be evenly distributed between them. The number 
of /„, n^, should not exceed the number of data Di we have, which is the total number of time slices, Nt. We choose 
the width of each Gaussian, Aw, to be half of the spacing between the Wj's. For a non-hermitian operator O, the 
strength function f{u)) is related to fi by 



^/.exp(-i(c.-c.,)VAV 



1 



(C9) 



while for an hermitian operator, f(—uj) = e f{<jj) in the canonical ensemble, and we choose 



1 



Aa;V27r 



(CIO) 



The response R{t) generated by /(w) is, for nonhermitian operators. 



(Cll) 



while for hermitian operators it is 



(C12) 



We choose the default model rrii, i = 1, n^^ to be a constant fixed by i?(0), which is related to the total strength. 
The data Dj are the values of R{t) for r = jA/3 measured from the Monte Carlo sampling. The error correlation 
function can also be measured in these calculations. 

To maximize the probability (C7), we have to know the dependence of Zg and Zq on a. Some simplification can 
be had by calculating Zs in the saddle point approximation. There is a saddle point in S{v?,m) at uf = rrii with the 
second derivatives \u^=mi — which leads to the approximate integral 



cfu exp{aS) ~ J (Fu exp{—2a'^^ul) = 



r/2 



(C13) 



The condition ^ 



then becomes 



1 f°° 

J-oo 



u S exp{aS - -x'^) 



1 dZs 
Zs da 



2a 



(C14) 



and the average image (fi) is given by 

1 r°° 

ifi) = ^ / d^n 

J-oo 



fiex.p{aS - -x^) 



(C15) 



We do these integrals by Monte Carlo sampling of Ui with exp(aS' — ^X^) as the weight function. The value a is 
determined by the self-consistent condition {S)a = When a is known, the average distribution {fi) = (u?) and 
the uncertainty Sfi = V {ff) — {fi)'^ is also given in the course of the Monte Carlo evaluation of (C15). 
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FIG. 1. Calculations in the grand canonical ensemble for protons only in the sd shell with monopole interaction (all six 
Ej=o matrix elements of the Wildenthal interaction), at (Np) = 3.17, (3 = 1- Shown are (H) and (J^) as functions of A/9 for 

three different decompositions: pure pairing decomposition, pure density decomposition, and a half-density and half-pairing 
decomposition. Solid diamonds at A/? = are the exact results obtained by direct diagonalization. 

FIG. 2. Zero-temperature calculations of ^''Mg with the schematic interaction (6.7). Note the relaxation of (H) and (J^) 
as /? increases. Hollow triangles are static path calculations in the pure density decomposition, solid diamonds are static path 
calculations by decomposing the pairing interaction into pair operators and the multipole interaction into density operators. 
Solid circles and hollow squares are both calculations in a pure density decomposition with A/3 = 0.125, using the Hartree 
solution and the maximal prolate state, respectively, as the trial wavefunction. The solid line segments indicate the exact 
ground state results. 

FIG. 3. Canonical ensemble calculations of ^°Ne with the schematic interaction (6.7) at A/3 = 0.25, 0.125 and 0.0625 and 
the exact results; {H) and { J^) are shown as functions of /3. These calculations were done in a pure density decomposition. 



FIG. 4. Similar to Fig. 3 for ^^Mg. 



FIG. 5. Finite temperature cranked calculations of '^'^Ne with the schematic interaction (6.7) in the canonical ensemble using 
a pure density decomposition. Here /3 = 1, with A/3 = 0.125, 0.0625 and 0.03125. The exact cranking curve is also shown. 



FIG. 6. Canonical ensemble calculations of the response functions for ^°Ne (/3 — 2.5) at discrete imaginary time using 
A/3 = 0.125,0.0625, in a pure density decomposition. The exact results are calculated in the ground state, (a), (c), (e) show 
the isoscalar quadrupole {Q = Qp + Qn), isovector quadrupole (Q„ = Qp — Qn) and the isovector angular angular momentum 
{Jv = Jp — J„) responses. The corresponding most probable strength functions recovered by the MaxEnt method are shown 
in (b), (d), (f) respectively. The exact strength functions calculated from ground state are plotted as discrete lines with the 
height indicating the integrated strength of the deltarfunctions. 

FIG. 7. Similar to Fig. 6. but for the single-particle pickup and stripping response, (a), (c), (e) show the imaginary 

time stripping response for the J = | orbital, and the pickup responses of the j = | and j = \ orbitals respectively. The 
corresponding most probable strength functions recovered by the MaxEnt method are shown in (b), (d), (f) respectively. The 
exact response and strength functions are calculated for the ground state. 

FIG. 8. Similar to Figs. 6,7 but showing the imaginary time pickup and stripping responses of j = | orbital. The response 
functions are in agreement with the exact curve for small t, but then abruptly follow the J = § response at larger r. 
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TABLE I. McuxEnt extraction of the moments of the strength functions corresponding to Figs. 6 and 7. The extrapolated 
(A/3 — » 0) total strength and first two moments are compared with the exact results for the ground state of ^"Ne. 







AO r\ IOC 


iXp — U.UDzo 


extrap 


exact 




total strength 


27.3 ± 0.2 


25.9 =b 0.1 


24.5 


25.1 


Q{t) ■ Q(0) 


(a;) 


2.33 ± 0.08 


2.77 ± 0.08 


3.22 


3.46 




(CJ ) 


O.Oy ± i.2 


10. ± i.2 


lz.9 


io.4 




total strength 


6.26 ± 0.03 


6.78 ± 0.02 


7.29 


6.96 


Qv{r)-Q^{Q) 


\U)/ 




7.77 ± 0.10 


8.31 


8.38 






59 9 ± 3 9 


66 6 ± 2 5 


73.4 


73.8 






16.3 ± 0.1 


1 6 05 + n 08 


15.8 


15.9 


Mr) ■ Jv{0) 


{uj) 


8.49 ± 0.25 


9.44 ±0.19 


10.39 


10.39 






89.8 ± 9.04 


107.7 ± 6.4 


125.6 


119.6 




total strength 


1.59 ±0.01 


1.62 ±0.07 


1.64 


1.59 


E^4/2mW«5/2m(0) 


(w) 


9.84 ± 0.12 


10.32 ±0.09 


10.80 


10.98 




98.0 ± 2 


107.5 ± 1.4 


117 


121 




total strength 


4.47 ±0.01 


4.42 ± 0.09 


4.37 


4.41 


E^OS/2m(T)4/2^(0) 




-3.15 ±0.02 


-3.00 ±0.02 


-2.86 


-2.81 




10.28 ± 0.05 


9.71 ± 0.04 


9.14 


10.08 




total strength 


1.702 ±0.004 


1.745 ± 0.003 


1.788 


1.773 




(co) 


-3.19 ±0.01 


-3.22 ± 0.02 


-3.25 


-3.16 




10.43 ± 0.06 


10.65 ±0.04 


10.87 


11.62 
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